Mon. Not. R. Astron. Soc. 000,[TJE5|(2005) Printed 5 February 2008 (MN I^TrX style file v2.2) 



A comparative study of disc-planet interaction 

M. de Val-Borro 1 *, R. G. Edgar 12 , P. Artymowicz 13 , P. Ciecielag 4 ' 5 , P. Cresswell 6 , 

G. D'Angelo 7 , E. J. Delgado-Donate 1 , G. Dirksen 8 , S. Fromang 6 9 , A. Gawryszczak 5 , 

H. Klahr 10 , W. Kley 8 , W. Lyra 11 , E Masset 12 ' 13 , G. Mellema 14 ' 1 , R. P. Nelson 6 , 
O S.-J. Paardekooper 14 , A. Peplinski 1 , A. Pierens 15 ' 6 , T. Plewa 16 , K. Rice 17 , C. Schafer 8 , 
(N R. Speith 8 

1 Stockholm University, AlbaNova University Center, SE-106 91, Stockholm, Sweden 

2 Dept. of Physics and Astronomy, University of Rochester, NY 14627, USA 

| 3 University of Toronto at Scarborough, 1265 Military Trail, Toronto, Ontario MIC 1A4, Canada 

, 4 University Observatory Munich, Scheinerstr. 1, D-81679 Munich, Germany 
' ■ 5 Nicolaus Copernicus Astronomical Centre, Bartycka 18, Warsaw, PL-00-716, Poland 
^ 6 Astronomy Unit, Queen Mary, University of London, Mile End Rd, London El 4NS, UK 

7 School of Physics, University of Exeter, Stocker Road, Exeter, EX4 4QL, UK 

8 Institute of Astronomy and Astrophysics Tubingen, Auf der Morgenstelle 10, D-72076 Tubingen, Germany 

9 DAMTP, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge, CB3 OWA, UK 



> 

1 10 Max-Planck-Institut fiir Astronomie, Konigstuhl 17, D-691 17 Heidelberg, Germany 



' 11 Department of Astronomy & Space Physics, Uppsala Astronomical Observatory, Box 515, 751 20, Sweden 
Y~ j 12 CEA, Service d'Astrophysique, Saclay, 91191 Gif-sur-Yvette Cedex, France 

13 IA-UNAM, Ciudad Universitaria, Apartado Postal 70-264, Mexico D.F. 04510, Mexico 



O 



^ 14 Leiden Observatory, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands 

| 15 Luth, Observatoire de Paris-Meudon, 92 195 Meudon Cedex, France 

^ r*| ! 16 ASC FLASH Center, University of Chicago, 5640 South Ellis, Chicago, IL 60637, USA 

Q j 17 Scottish Universities Physics Alliance, Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh, EH9 3HJ, UK 

6 ■ 



5 February 2008 



x 



ABSTRACT 

We perform numerical simulations of a disc-planet system using various grid-based and 
smoothed particle hydrodynamics (SPH) codes. The tests are run for a simple setup where 
Jupiter and Neptune mass planets on a circular orbit open a gap in a protoplanetary disc dur- 
ing a few hundred orbital periods. We compare the surface density contours, potential vorticity 
and smoothed radial profiles at several times. The disc mass and gravitational torque time evo- 
lution are analyzed with high temporal resolution. There is overall consistency between the 
codes. The density profiles agree within about 5% for the Eulerian simulations while the SPH 
results predict the correct shape of the gap although have less resolution in the low density 
regions and weaker planetary wakes. The disc masses after 200 orbital periods agree within 
10%. The spread is larger in the tidal torques acting on the planet which agree within a factor 
2 at the end of the simulation. In the Neptune case the dispersion in the torques is greater than 
for Jupiter, possibly owing to the contribution from the not completely cleared region close to 
the planet. 

Key words: Physical data and processes: accretion, accretion discs - hydrodynamics. Solar 
system: planets and satellites: general. 



1 INTRODUCTION which can be avoided, due to the central part that gas plays in the 

cosmos. 

Hydrodynamics is a difficult subject, which has caused many prob- 
lems for many distinguished physicists. However, it is not a topic 
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The basic equations of hydrodynamics are the Navier-Stokes 
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equations, and have been known for almost two centuries: 

^ + V-(pv) = (1) 

^ + (v-V)v = --V/?-V<£ + V-T (2) 
dt v J p 

where p is the density, v the velocity of the fluid, p the pressure, 
the gra vitational potential and T is the fu ll viscous stress tensor 
(see e.g. Mihalas and Weibel Mihal alll984l) . The first equation de- 
scribes the conservation of mass and the second, conservation of 
momentum. An equation of state closes the system of equations, 
and additional terms may be added as required. Despite their com- 
paratively simple form, the Navier-Stokes equations have proved 
remarkably stubborn to mathematical analysis. 

The problem lies in the v • Vv terms (the so-called 'advection' 
terms). These arise because the equations describe a fluid moving 
past a fixed point in space (the Eulerian point of view). The ad- 
vection terms make the equations non-linear (since they are effec- 
tively proportional to v 2 ), rendering many mathematical techniques 
useless. Indeed, no one has yet proven that solutions to the Navier- 
Stokes equations are unique. This is in sharp contrast to many other 
important equations in physics. For example, the Poisson equation 
is linear, and has unique solutions. This opens up many avenues for 
obtaining solutions — the Method of Images being a well known 
example. As a result of this non-linearity, theoretical investigation 
of fluids has to be restricted to highly idealised flows. 

To make progress then, we are forced to turn to computers. 
Numerical algorithms for solving complex equations have been 
studied for centuries, and computers are ideal for implementing 
these. Unfortunately, computers are tricky beasts, with a habit of 
doing precisely what you told them to, just when you least ex- 
pected it. The 'obvious' way of computing a numerical solution 
may well be unstable (this is particularly true of the Navier-Stokes 
equations), and an implementation of a stable method may well 
contain bugs. Floating point numbers have finite accuracy, and var- 
ious subtleties arise when codes approach this limit. In particular, 
arithmetic ceases to be distributive and tests for equality cease to 
be reliable. Different architectures, operating systems and compil- 
ers all add to the mix. For this reason, work performed on com- 
puters is better described as a 'numerical experiment' rather than a 
'simulation.' 

Fortunately, there is no need to be overly pessimistic about 
the situation. For example, although the minutiae of floating point 
numbers and the vagaries of different compilers can be trouble- 
some, these should not give problems in the majority of cases. Un- 
less agreement to the last bit is required, there should be no sig- 
nificant difference in results obtained with different architectures 
and compilers. Instead, it suffices to focus on differences between 
algorithms for solving the Navier-Stokes equations. Of course, all 
hydrodynamics codes are carefully tested against simple problems 
(such as shock- tubes). It is on more complex problems that differ- 
ences and difficulties can be exposed. 

Within the context of the EU-RTN "The Origin of Planetary 
Systems," 1 we have conducted a comparison of hydrodynamics 
codes, which we present in this paper. The problem we selected was 
that of a planet in a fixed circular orbit in a circumstellar disc. This 
has the virtue of simplicity, while still retaining sufficient complex- 
ity to allow us to see meaningful differences between the various 
algorithms. We ran the test problem on 17 independent codes. 
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A comparison of several numerical methods o n the pro blem 
of a p lanet embedded in a disc was performed bv iBrvden et alJ 
( 1999) using SPH, van Leer and Godunov methods with different 
equations of state. In particular, they studied the accretion onto the 
planet after it had cleared a gap. Other examples of comparisons in 
different fields to verify algorithms and implementations published 
during the last few years include the Santa Barbara cluster project 
( Frenk et al. 1999), the non-LTE radiative transfer code comparison 
(Ivan Zadelhoff e t al. 2002), the Ra yleigh-Taylor i nstability study 
by the Alpha-Group collaboration (Di monte et alJl2004|) and the 
comp arison of models of photoionization regions iPeauignot ^etal] 
l200lh . 

The aim of this project is to test the reliability of present nu- 
merical computations of disc-planet interaction with a quantitative 
comparison and generate a benchmark for future simulations. In 
Section|2| we briefly describe the interaction between a planet and 
a protoplanetary disc and outline the motivation for this study. The 
initial setup and boundary conditions of the problem are described 
in Section [3] In Section |4| the numerical methods used in the com- 
parison are described. The results are shown in Section [5] We dis- 
cuss the results in Section[6] and in Appendix IaI we summarize our 
experience with this project that could be useful for future compar- 
isons. 



2 DISC-PLANET INTERACTION 

Over 160 extrasolar planetary systems have been discovered 
by radial velocity measurements during the l ast years (e.g. 
Mav or and OuelozN l995: Ma rcv and ButleJ h996). Giant planets 
have been found in very close orbits around the central star with 
orbital periods of a few days and almost circular orbits, the so- 
called "Hot Jupiters". Planets orbiting at larger distances from the 
star show a broad eccentricity distribution reaching roughly e = 0.9 
(for recent review s of th e properties of the observed systems see 
Marc v et alJl2003l 12005). The origin of the differences with the 
planets in the solar system is not well understood, although various 
explanations have been proposed. The standard models explain gi- 
ant planet formation either through planetesimal acc umulation fol- 
lowed by rapid gas accretion onto the planet core ( Pollac k et al J 
1996) or gravitational instabilities in the disc (see e.g. lBossl 1998. 
2001). In both cases the planets are likely to have formed at larger 
distances from the central star than observed. 

Orbital migration due to gravitational interaction between the 
planet and the gaseous disc is a possible mechanism to bring 
planets to a close orbit. The tidal interaction between a planet 
and a gaseous disc wa s studied before the discovery of ex traso- 
lar planetary system s by Goldreich and Tremaind o 1 979 . 1980) and 
iLin and Papaloizoul 1 19791 Il986allbl) . In the linear approximation 
the planet excites waves at the Lindblad resonances that deposit 
angular momentum in the disc. The flux of angular momentum has 
different signs in the inner and outer discs causing the orbital mi- 
grati on of the plan et. 

IWard dl997h proposed that two different types of planetary 
drift exist. Type I migration occurs when the planet mass is small 
and migrates relative to the disc with a rate proportional to its mass 
and the surface density of the disc. This migration is quite fast and 
the orbital decay timescale of the order of 10 5 years is comparable 
to the formation timescale of a giant planet by planetesimal accu- 
mulation. In type II migration the planet is massive enough to open 
a gap in the disc. The planet is then locked to the viscous evolu- 
tion of the disc and its migration rate will be determined by the 
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strength of the viscosity. The estimated timescale for type II migra- 
tion is one or two orders of magnitude larger than the type I mi- 
gration timescale for the same planetary mass. Type II migration is 
believed t o be responsible for the presence of p lanets at short orbital 
distances dTrilling et alJl200lludrv et all20oA . Nu merical simu- 
lations of planet migration i n a viscous disc (see e.g. iNelson et alJ 
2000; D'Angel o et alJEooih confirm the inward migration of the 
planet on the viscous time- scale predicted by linear theory for both 
accreting and non- accreting planets. 

The nonlinear interaction between disc and planet cannot be 
fully described analytically or reproduced in laboratory experi- 
ments. Therefore, multidimensional hydrodynamical simulations 
of protoplanetary discs with embedded planets during many or- 
bital periods are necessary to understand the formation and evo- 
lution of extrasolar planetary systems. However, some differences 
are found in the simulations depending on the numerical algo- 
rithm employed. The spiral waves generated around the planet may 
be stationary in the c o-rotating frame (e.g. ZEUS -based results 
of iLubow et alJll99^ . Other higher-order hydrodynamical codes 
show time variability of the flow in the spiral arms propagating 
along the shock. The quasi-periodic disturbances in the shocks have 
important implications for the formation and evolution of vortices 
along the edges of the gap opened by the planet. In some simula- 
tions wavy structures and vortices are observed at the edge of the 
gap opened by the plan et which interact with the shocks (see e.g. 
Nelson and Benz l2003h . In this paper, we have used different algo- 
rithms presently in use in the astrophysical community to study the 
planet-disc system in a simple but meaningful case. 



3 SETUP DESCRIPTION 

We examined the gap opening by a giant planet in an infinitesimally 
thin disk with a constant surface density. The numerical setup was 
defined in the web 2 where interested modelers were invited to par- 
ticipate in the comparison. 

The planet's gravitational potential was given by the formula 



vV 2 + £ 2 



(3) 



where r is the distance from the planet and £ is the gravity soften- 
ing. The simulations were run with two different softening coeffi- 
cients: 



£l = 0.2r L 



(4) 



where = (/i/3) 1 / 3 is the size of the Roche Lobe of the planet, 
and the larger value 



£ 2 = 0.6#j 



P' 



(5) 



with Hp the disc scale height at the planet location. The second 
softening was mainly introduced to mimic the torque cut-off due 
to the effect of the disc vertical distribution. The results discussed 
in this paper concern mostly the calculations that use the larger 
softening. In our simulations, the self-gravity, energy transfer and 
magnetic fields in the gaseous disc were not considered. 

The mass relation between the planet and the star was chosen 
so that the reduced mass has the values fi = M p /(M* +M p ) = 10 -3 
and 10 -4 , corresponding to roughly Jupiter and Neptune masses 
when the star mass M* = M . The planet was kept in a circular 
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orbit at approximately semi-major axis a = 1 ignoring the effect of 
the gravitational torques on the planet. The position of the planet 
in the cell with respect to the cell's corner is given in Table 1. 
The computations were performed in the radial domain [0.4a, 2.5a] 
to study the influence of the planet in a sufficiently large fraction 
of the disc. In the cylindrical grid codes, the number of cells in 
the radial and azimuthal directions were n Y xn^ = (128,384) with 
uniform spacing in both dimensions. Therefore, the cells around 
the planet position were approximately square. Several tests were 
done with different schemes at resolution n Y xn^ = (256,768) and 
nrXrifj) = (512, 1536) to check the convergence of the results. The 
polar coordinates schemes used a corotating reference frame. The 
centre of the frame was not specified in the problem description 
and codes with frames centred in the centre of mass (CM) and cen- 
tral star were used. The star position was fixed at r, = (0, 0) and 
the planet at r, (j) = (1,0) in co-rotating coordinates, where the az- 
imuthal range was [— tt, tt]. The Cartesian schemes FLASH-AP 
and Pencil were run on a uniform non-rotating grid at resolution 
n x xn y = (320, 320), and n x xn y = (640, 640). The computational 
domain was [—2.6a, 2.6a] x [—2.6a, 2.6a]. The unit of time used in 
the simulations was the orbital period at a = 1 which is defined as 



P v = 2n 



.1/2 



G(M*+M p ) 



:27T, 



(6) 



where G = 1 and M* +M p = 1. The angular frequency of the planet 
was £l D = 1 in our units. 



3.1 Initial conditions 

The modelled disc was 2-dimensional so that the vertically inte- 
grated quantities were solved. The initial surface density was con- 
stant and given by the expression 



: 0.002 - 



9 



(7) 



where a is the semi-major axis of the planet. We assume that the 
heat generated by viscous dissipation and tidal forces in the disc is 
radiated away, so the disc remains geometrically thin. The initial 
angular velocity was fixed to the local Keplerian frequency at the 
given radial position and the radial velocity was zero initially. 

We used the standard sound speed profile of a slightly flaring 
solar nebula H/R = c s /vk = 0.05, where H is the disc scale height, 
R the distance from the centre of the star, vk the local Keplerian 
velocity, and c s the isothermal sound speed defined as 



dp 
dJL' 



(8) 



which has a dependence on radius c s oc r -1 / 2 . This corresponds 
to a locally isothermal equation of state with a profile T(r) oc r~ l 
maintained through the simulation. The disc height at the planet 
location remains constant during the opening of the gap. 

The planet mass was gradually increased during the first 5 or- 
bital periods using the expression 



M(t) 



Kt 

10ft 



(9) 



to avoid the the appearance of strong shocks seen when the planet 
is introduced instantaneously. The gas accretion from the disc onto 
the planet was ignored. This situation can be realistic in the case 
when the planet's atmosphere fills the Roche lobe and no further 
accretion is allowed. 
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The problem was originally proposed to be run with no artifi- 
cial viscosity or as low as possible as allowed by the code. Some of 
the used codes include artificial viscosity to smooth out the shock 
fronts and prevent unphysical results as described in Section|4| We 
performed simulations for each planet mass including a physical 
viscosity that genera tes a stress tensor with a turbulent vi scosity co- 
efficient v (see e.g. lLandau and Lifsh itz 1959; Kiev 1999). The val- 
ues of the kinematic viscosity used in our simulations were v = 
and 10 -5 in units where a = 1 and G(M* +M p ) = 1. The simula- 
tions were run typically during several hundred orbital periods for 
each of the planet masses and viscosity coefficients. 



3.2 Boundary conditions 

To completely define the problem, we describe the implemented 
boundary conditions. The disc was considered as an isolated system 
with no mass inflow. We used solid boundary conditions with wave 
killing zones next to the boundaries to reduce wave reflection in the 
cylindrical coordinates codes. In the polar coordinates schemes, the 
damping regions were implemented in the radial ranges [0.4a, 0.5a] 
and [2. la, 2.5a], where the following equation was solved after 
each timestep: 



dx 
dt 



-XQ 



R(r) 



(10) 



where x represents the surface density and velocity components, t 
is the orbital period at the corresponding boundary and R(r) is a 
parabolic function which is one at the domain boundary and zero at 
the interior boundary of the wave killing zones. This wave damping 
condition does not conserve mass, but the mass loss is very small 
as shown below. 

The grid-based codes in Cartesian coordinates implemented 
the same wave killing condition as the polar codes in the ring 
[2. la, 2.5a]. Tests were done including the damping condition in 
the region [0.4a, 0.5a] although this is not necessary since there 
is no inner solid boundary. There was free outflow in the x and y 
boundaries in the Cartesian implementations. 

Note that the SPH codes implement different boundary condi- 
tions using rings of virtual particles as described in section[4] 



3.3 Output data 

2-D snapshots of the density and velocity components were output 
at 2, 5, 10, 20, 50, 100 and 200 orbital periods for grid codes, al- 
though in some cases the simulations were run up to 500 periods. 
All the physical quantities were given at the cell centres. 

In the case of SPH codes the output quantities at the previ- 
ous times were particle positions, velocity components, smoothing 
length and mass. The particles were projected to a 2-dimensional 
cylindrical grid with the resolution n v x n§ = (128,384) to com- 
pare directly with the lower resolution results from the Eulerian 
grid codes. The associated kernel for each particle used internally 
by our codes was the th i rd o rder spline function introduced by 
Monagha n and Lattanziol ll985h with a multiplicative coefficient 
corresponding to a 2-dimensional simulation. The density at a given 
point was calculated by interpolation with the spline kernel using 
the expression 



<p(n)> ^ 



N 

: £mjW( 

7=1 



(11) 



and |rj — rj | is the distance from the cell centre to the given parti- 
cle. The smoothing length hj has different values for each particle. 
In a similar manner, the velocity components were interpolated to 
the grid with the kernel function and normalized with respect to 
the integrated kernel. The resolution element of the SPH models 
is given by the smoothing length of the particles. For the number 
of particles used in the SPH calculations, the effective resolution 
is similar to the number of cells in the hydro models at the afore- 
mentioned resolution, were the particles distributed in an equiva- 
lent spatial domain. SPHTREE uses a smaller smoothing length 
than PARAS PH and therefore should have a slightly better spatial 
resolution in our calculations. 

The azimuthally averaged density was obtained as 



1 f 2 * 
2-Jo Ld * 



(12) 



Slices of the surface density were taken at the planet position and 
Lagrangian points in the radial and azimuthal directions. 

We calculated the vortensity or potential vorticity, defined as 
the ratio of vorticity and surface density 



(Vxv) 



(13) 



In the frame rotating with the planet, the vortensity is given by the 
expression (V x v + 2£2 p )/Z, where £2 p is the orbital frequency of 
the planet. 

The gaseous disc interacts gravitationally with the planet 
by means of the torques generated by the spiral arms ( see e.g. 
Goldreich and Tremaine 1979; Papaloi zou and Linlll984f) . Every 
few timesteps the contributions from the inner disc excluding the 
Hill sphere, outer disc excluding the Hill sphere and the torque 
from the material between 0.5 and 1 Hill radius to the torque are 
recorded. The disc mass interior and exterior to the planet orbit 
was also obtained with the same output frequency. 

The torques were calculated in units where a — 1, P = 2% and 
M* = 1 — ji integrating over the corresponding region. In the case 
of a 2-dimensional disc the torque has only a vertical component 
which is given by 



GMr 



Er p x 



( r 2 + £ 2)3/2 



rdrd(j) 



(14) 



where raj is the mass of the particle, W(r,hj) is the spline kernel 



where L is the surface density, r p is the planet position, r e is the 
distance between the planet and the fluid element. 

We performed Fourier analysis of the torque data to under- 
stand the cause of the observed var iability. We used a Welch win- 
dowing function JPress etalJll992h to smooth the deviation be- 
tween the initial and final amplitudes in the time series. 



4 DESCRIPTION OF THE CODES 

We will now discuss the codes used in the comparison. Even within 
the restricted field of astrophysical fluids, there are many different 
algorithms for computing flows. There are then different imple- 
mentations of the same algorithm. We will therefore start with a 
discussion of the general principles of various types of codes pre- 
sented in this paper, and then go on to detail particulars of each 
code used. This is not meant to be a general review of all the types 
of codes used to conduct numerical experiments in astrophysics. 
For more detailed information, the read er should refer to any of 
the plethora of books on the subject (e.g. 
lLeVeauel200l . 

The parameters of each code are given in Table 1, including 



Table 1. Summary of the parameters used in all codes. Column 2: name of the users of the code; col. 3: reference to detailed code description; col. 4: numerical method (upwind, high-order 
finite-difference, shock-capturing or SPH); col. 5: Courant number used; col. 6: type of artificial viscosity used (none, von Neumann-Richtmyer, tensor or scalar); col. 7: reference frame of 
hydrodynamic codes (corotating or inertial); col. 8: centre of reference frame (centre of mass of star-planet system or primary star); col. 9: position of the planet in the cell (centre, corner, 
arbitrary or coordinates with respect to the cell's corner in units of the radial and azimuthal steps dr and d(j)); col. 10: processor; col. 11: approximate execution time on a single processor in 
hours per 100 orbits. 



Name 


Users 


References 


Method 


Courant 


Art. vise. 


Frame 


Centre 


Planet 


Processor 


Tcpu 


NlRVANA-GDA 


G. D'Angelo 


1 


Upwind 


0.5 


none 


corot. 


CM 


corner 


Power5 1.65 GHz 


4 


NlRVANA-GD 


G. Dirksen 


1 


Upwind 


0.67 


none 


corot. 


CM 


(arb.,0) 


P4 2.8 GHz 


6 


NlRVANA-PC 


P. Cresswell 


1 


Upwind 


0.5 


none 


corot. 


primary 


(0.29,0) 


P4 2.4 GHz 


11 


RH2D 


W. Kley 


2 


Upwind 


0.75 


1.0 (bulk) 


corot. 


primary 


(arb.,0) 


P4 3 GHz 


3.4 


GLOBAL 


S. Fromang 


3 


Upwind 


0.5 


1.0 


corot. 


primary 


(arb.,0) 


Xeon 3 GHz 


16 


FARGO 


F. Mas set 


4,5 


Upwind 


0.5 


2.0 (vN-R) 


corot. 


primary 


(0.57,0.5) 


P4 2.8 GHz 


1.25 


GENESIS 


A. Pierens 


4,5 


Upwind 


0.5 


1.0 (tensor) 


corot. 


primary 


centre 


P4 2.8 GHz 


1 


TRAMP-vanLeer 


H. Klahr & W. Kley 


6 


Upwind 


0.4 


1.1 (vN-R) 


corot. 


CM 


arb. 


Opteron 2 GHz 


16 


Pencil 


W. Lyra 


7 


High-order fin.-diff. 


0.4 


1.0 (bulk) 


inertial 


CM 


arb. 


P4 2.4 GHz 


36 


AMRA 


P. Ciecielag & T. Plewa 


8 


Shock-capturing 


0.8 


none 


corot. 


primary 


(0.57,0) 


Xeon 3 GHz 


21 


FLASH-AG 


A. Gawryszczak 


9 


Shock-capturing 


0.8 


none 


corot. 


primary 


(0.57,0.5) 


Athlon 2 GHz 


42 


FLASH-AP 


A. Peplinski 


9 


Shock-capturing 


0.7 


none 


inertial 


CM 


arb. 


Athlon 1.8 GHz 




TRAMP-PPM 


H. Klahr 


10,11 


Shock-capturing 


0.8 


none 


corot. 


primary 


arb. 


Opteron 2 GHz 


28 


Rodeo 


S-J. Paardekooper & G. Mellema 


12 


Shock-capturing 


0.8 


none 


corot. 


primary 


arb. 


Athlon 1.7 GHz 


25 


JUPITER 


F. Mas set 


13 


Shock-capturing 


0.7 


none 


corot. 


primary 


(0.57,0.5) 


P4 2.8 GHz 


15.3 


SPHTREE 


K. Rice 


14 


SPH 


none 


bulk + shear 




CM 


arb. 


Opteron 1.8 GHz 


10 


PARASPH 


C. Schafer & R. Speith 


15 


SPH 


none 


bulk 




CM 


arb. 


Opteron 2 Ghz 


250 



R eferences: l: | Zi egler and YorkeNl997T): 2: 1 
8: IPlewa and Mulled fcOOllh 9: iFrvxell et al J 
14: lBenzU l990h 15: Schaf er etalJfcOOA . 



\): 3:lHawlevandSt. 
10: Blondin and Lufki 



w , 4: iMasseJ fcOOOah: 5: iMasseJ bOOOhh: 6: iKlahr et al J J 1 999h: 7: [Brandenburg and Dobler (2002) 
; 11: IColella and Woodward il984h : 12: IPaardekooner and Mellem3 J2006h : 13: IPember etaTI 1 1995 s ) 
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Table 2. Codes that were run for the lower resolution runs of the setups 
defined in Section [3] The grid size for the Eulerian codes was n r x n§ = 
(128,384) for the cylindrical grid codes and n x x n y = (320,320) for 
FLASH- AP and Pencil. The number of particles was 250000 in SPH- 
TREE and 300 000 in PARASPH. 





Jupiter 


Jupiter 


Neptune 


Neptune 


Codes 


inviscid 


viscous 


inviscid 


viscous 


Nirvana- GD A 


X 


X 


X 


X 


Nirvana- GD 


X 


X 


X 


X 


Nirvana-PC 


X 


X 


X 


X 


RH2D 


X 


X 


X 


X 


GLOBAL 


X 


X 


X 


X 


FARGO 


X 


X 


X 


X 


GENESIS 


X 


X 


X 


X 


TRAMP-vanLeer 


X 




X 




Pencil 




X 




X 


AMRA 


X 


X 


X 


X 


FLASH-AG 


X 


X 


X 


X 


FLASH-AP 


X 


X 


X 


X 


TRAMP-PPM 


X 




X 




Rodeo 


X 


X 


X 


X 


JUPITER 


X 


X 


X 


X 


SPHTREE 




X 




X 


PARASPH 


X 


X 


X 


X 



Table 3. Codes that were run at resolution n r x n$ = (256,768) and equiv- 
alent resolutions for the Cartesian grid and SPH schemes. 





Jupiter 


Jupiter 


Neptune 


Neptune 


Codes 


inviscid 


viscous 


inviscid 


viscous 


Nirvana- GD 




x 






Nirvana-PC 


x 


x 


X 




AMRA 


x 


X 






FLASH-AP 


x 


X 


X 


X 


PARASPH 




X 







references in which the algorithms are described in detail. Table 
shows which codes were run for the low resolution defined tests. 
In Tables [5] and |4] we show the schemes that were run at higher 
resolution. 

4.1 Grid Based Codes 

As the name implies, grid based codes cover the computational vol- 
ume with a set of grid points at which the various flow variables 
(velocity, density etc.) are computed. The mesh geometry (conven- 
tionally orthogonal, although this is not absolutely required) can be 
chosen to reflect the underlying symmetry of the problem. This of- 
ten leads to a reduction in the number of grid cells required for 
a particular problem, and a corresponding cut in computational 
time. The codes used in our problem use a reference frame cen- 

Table 4. Codes that were run at resolution n x xn§ = (512, 1536). 





Jupiter 


Jupiter 


Neptune 


Neptune 


Codes 


inviscid 


viscous 


inviscid 


viscous 


Nirvana- GD 




X 






RH2D 




X 






FARGO 


X 


X 


X 


X 



tred in the centre of mass or primary as indicated in column 8 of 
Table 1 . All the simulations centred on the primary include the indi- 
rect terms in the potential. For astrophysical (compressible flow at 
high Reynolds number) flows, two different approaches to solving 
the fluid equations are generally used. However, before we describe 
these, some general points should be noted. 

The most important of these is the Courant-Friedrichs-Lewy 
(CFL) condition. Simply stated, informatio n must not t ravel more 
than one grid cell per timestep (see, e.g. IPress et all dl992h for 
a mathematical derivation). In a hydrodynamics code, this trans- 
lates into a restriction on the timestep, based on ve locity and sound 
speed (some authors, e.g. Edgar and Clarke (2004), have also added 
an acceleration condition when appropriate). Violation of the CFL 
condition leads to unphysical effects, as causality is violated. When 
we refer to the "Courant number" in the descriptions below, we are 
describing an extra safety factor, beyond the formal CFL condition 
itself. Note, however, that the CFL condition only applies to time- 
explicit codes. Implicit solvers are not restricted by it, but no results 
based on such a code were submitted to us. 

Next is the extension to multi-dimensions. Most algorithms 
for solving the equations of hydrodynamics have been developed 
for one dimensional flow. The conventional method for using a one 
dime nsional algorithm in multiple dimensions is Strang splitting 
lStranJl968l) : solve the ID equations along each row of cells (the 
x\ direction), then solve along each column (thex2 direction), using 
the updated values from the x\ sweep. Formally, the x\ step should 
be split in two as ^x\ —> X2 —> \x\ , but most codes do a full step in 
each direction and alternate which is done first (this is sometimes 
called "Godunov splitting"). The Strang approach makes orthogo- 
nal co-ordinates highly desirable. To minimise the truncation errors 
this approach produces, the grid cells must be kept locally square. 

Most codes presented here use a rotating polar grid. For these, 
there is an extra subtlety: the treatment of the Coriolis force. As is 
conventional in fluid dynamics, the simple and obvious way to in- 
clude this (as an extra force) leads to incorrect angula r momentum 
transport. Instead, the angular momentum approach of lKlevH l998) 
must be used. On reflection, this is unsurprising: the Coriolis force 
simply enforces the conservation of angular momentum in a rotat- 
ing frame. 

Although not relevant to the comparison problem itself, many 
of the codes here can make use of refined meshes. High resolution 
is always desirable, but computationally expensive. To concentrate 
grid cells where they are needed, patches of the grid may be cal- 
culated at higher resolution, and the results communicated back to 
the coarser parent grid. Patches can themselves be patched, giving 
the potential for extremely high resolution. If the patches are de- 
termined at the start of a calculation, such a code is said to be of 
the 'nested grid' type. However, some codes can dynamically add 
and remove patches. This is known as adaptive mesh refinement 
(AMR). For this comparison, we have chosen not to use refined 
meshes. This is in the interests of simplicity, since there are a vari- 
ety of algorithms for performing the refinement, and we are already 
comparing a large number of codes. However, we would encourage 
other workers in the field to compare refinement methods. 

4.1.1 U pwind Methods 

The upwind codes used in this comparison work by discretising 
the appropriate version of the Navier-Stokes equations, and solving 
that. These codes use the technique of operator splitting, and some 
operators are discretised in a finite difference manner, while others 
are solved with a finite volume method. For this reason, codes simi- 
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lar to those we shall now discuss are sometimes refered to as 'finite 
difference/volume' schemes, or even just 'finite difference.' We es- 
chew this epithet, since almost any grid code could be described as 
'finite difference' at some level. 

In a typical operator split scheme, each timestep is split into 
two phases. During the source step, the velocity is updated using 
the source terms in the Navier-Stokes equations (pressure gradi- 
ents, gravity etc.). In the transport step, these velocities are then 
used to advect (the v • Vv terms) the other quantities. This is usu- 
ally done conservatively using the integral form of the equations 
(integrated over a volume - hence the name). During the advection 
step, second order 'upwinding' is used (interpolation based on ve- 
locities), to ensure that shocks remain sharp. Some sort of artificial 
viscosity is generally required to stop post- shock oscillations mak- 
ing the code unstable. 

These codes usually use a staggered mesh, to improve the 
order of their differencing. Scalar variables (such as density) are 
stored at zone centres, while vector quantities are stored at the faces 
(e.g. vi is stored at the centre of the x\ face). 

Codes like these are often des cribed as being 'ZEUS-l ike' - 
a reference to the ZEUS code of Stone and Norman ( 1992). Al- 
though that paper provides an excellent description of the methods 
used, the epithet 'ZEUS-like' does not generally mean "derived 
from ZEUS." Rather, they are based on the same or similar algo- 
rithms, and ZEUS happens to be the best known implementation 
of these. 

4.1.1.1 The Nirvana code In this comparison, three sets of 
results were submitted w hich made use of the NIRVANA code of 
Ziesler and Yorke (1997). All of the following codes are based on 
the original version of NIRVANA, which was not publically re- 
leased. Each of the codes was enhanced from the original code 
base by different groups over a number of years. Hence, variations 
between the NIRVANA codes highlight how even the same basic 
algorithm can vary. Different Courant numbers were also used - 
Nirvana- GD used 2/3, while Nirvana- GDA and Nirvana-PC 
used 1 /2. 

4.1.1.2 The RH2D code The RH2D code is two-dimensional 
mixed explicit/implicit 2nd order upwind algorithm on a staggered 
grid. The adve ction algo r ithm is based on the monotonic trans- 
port scheme bv lvan Leea Jl977h . The RH2D code can treat radi- 
ation transport in the flux-limited diffusion approximation, and in- 
cludes the full tensor viscosity with dissipation. In contrast to some 
other codes the velocity variables that are evolved in RH2D are 
radial v and angular velocity £L Both radiation and viscosity can 
be solved impli citly to avoid possible time-step limitations. We re- 
fer the reader to lKlevI Jl989h for a full description of the code. For 
the purpose of the present calculations the radiation module was 
replaced by a locally isothermal equation of state. The viscosity 
was solved explicitly. The formulation of the equations, in particu- 
lar the treatment of the physical and artificial viscosity in the stress 
tensor components, has been d escribed with respect to the embed- 
ded planet problem in detail by Kiev ( 1999). 

4.1.1.3 The GL O BAL code The GLOBAL code 
jHawlev and Stond Il995h is derived from ZEUS 
( Ston e and NormarJ 1 19921) . The Courant number was 0.5, and 
an artificial viscosity co-efficient of 1.0 was required to stabilise 
wave propagation in the disc. 



4.1.1.4 The FARGO code FARGO is a simple 2D polar mesh 
code dedicated to disc planet interactions 3 . It is based upon a stan- 
dard, ZEUS-like hydrodynamic solver, but owes its name to the 
FARGO algorithm upon which the azimuthal advection is based 
(Masset 2000a b). This algorithm avoids the restrictive timestep 
typically imposed by the rapidly rotating inner regions of the disc, 
by permitting each annulus of cells to rotate at its local Keple- 
rian velocity and stitching the results together again at the end of 
the timestep. The use of the Fargo algorithm typically lifts the 
timestep by an order of magnitude, and therefore speeds up the 
calculation accordingly. The mesh centre lies at the primary, so in- 
direct terms coming from the planets and the disk are included in 
the potential calculation. The Courant number was 0.5, and a sec- 
ond order artificial viscosity of = 2 (cf equations 33 and 34 of 
IStone and^N orman) was used. 

The standard boundary conditions prescribed in the test prob- 
lem were used. In addition, the dependence of the results on the 
damping condition was tested using a slightly different boundary 
were a transmitted wave boundary condition was used. The pitch 
angle of the wake at the inner and outer boundary was valuated us- 
ing the WKB approximation. The content of the border ring was 
then copied into the ghost ring, properly azimuthally shifted by the 
amount dictated by the pitch angle. This technique is very efficient 
at removing any reflected wave and yields similar results to the 
standard boundary condition defined in Section I3T2I 

4.1.1.5 The GENESIS code GENESIS is a 2D code which 
solves the fluid equations using a upwind method with a time- 
explicit, operator- splitting procedure. The FARGO algorithm (see 
description above) is applied to avoid the timestep limitation at the 
inner edge of the disc. Because of this, the code does not alternate 
radial and azimuthal integrations. Artificial viscosity is handled by 
using a bulk viscosity in the viscous stress tensor dKlevI 19991) . 

4.1.1.6 The TRAMP-vanLeer code This is a 3D version of 
RH2D (see above) with the same second order van-L eer scheme 
(simila r to that used in the ZEUS and NIRVANA codes \ lKlahr et alJ 
( 1999) provide a description. The fact that the code is intrinsically 
3D explains why it performs two times slower than the pure 2d ver- 
sion RH2D 4 . We use a moderate value of 1 . 1 for the von Neumann- 
Richtmyer type viscosity. The implementation works in the corotat- 
ing frame where the centre of mass is the centre of the coordinate 
system. Hence no extra acceleration terms are necessary. 

4.1.2 High- Order Finite -Difference Methods 

4.1.2.1 The Pencil code Pencil is a non-conservative finite- 
difference code that uses sixth order centred spatial derivatives and 
a third order Runge-Kutta time- stepping scheme, being primarily 
designed to deal with compressible turbulent magnetohydrodynam- 
ical flows 5 . Being high-order, Pencil needs viscosity and diffusiv- 
ity terms in order to stabilize the numerical scheme. For this reason, 
we could not perform inviscid runs. 

The code is intrinsically 3D and Cartesian, structured in a 
cache-efficient way. The domain is tiled in the y and z direction for 
parallelization, with the original 3D quantities being split into ID 

3 FARGO is available at http://www.maths.qmul.ac.uk/~masset/fargo/ 

4 The remaining factor of two comes from the roughly two times smaller 
Courant number in TRAMP-vanLeer 

5 Pencil is available at http://www.nordita.dk/software/pencil-code/ 
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arrays - pencils - in the x direction, hence the name of the code. The 
equations are solved along these pencils in the x direction, which 
leads to the convenient side-effect that auxiliary and derived vari- 
ables use very little memory as they are only ever defined on one 
pencil. By calculating an entire timestep in the x-direction along 
the box, Pencil can achieve a speed-up of ~60% on typical Linux 
architectures. 

This is the first time Pencil has been applied to the embedded 
planet problem. 



4.1.3 Shock- Capturing Methods 

The other scheme for grid- based astrophysic al fluid flows in com- 
mon use is that proposed bv lGodunmi Jl959h . Such schemes make 
use of the fact that there is an analytic solution to the ID shock 
tube problem: the so-called Riemann problem. Godunov's original 
scheme treated each cell as piecewise-constant (i.e. variables such 
as density were assumed to be c onstant throughout the cell), g iving 
a sharp shock at each interface. IColella and Woodw ard ( 1984[) im- 
proved Godunov's method by using parabolic interpolation, giving 
the 'piecewise parabolic method' (PPM) which is the most com- 
mon implementation in use today. Implementations of PPM can 
be in Eulerian or Lagrangian form. For the purposes of the inter- 
polation, all values are stored at the cell centres (cf the staggered 
grids mentioned above). Shock-capturing codes include the pres- 
sure gradient in the basic solver. Since solving the full Riemann 
problem is computationally expensive, many codes use an approx- 
imate solver. Furthermore, to deal with strictly isothermal flows, a 
special isothermal Riemann solver must be written, since the con- 
ventional one involves 7 — 1 denominators. 

Shock-capturing schemes do not usually require any artifi- 
cial viscosity to ensure stability (sometimes authors will include 
a small artificial viscosity to prevent post-shock oscillations, but 
these oscillations do not usually threaten the stability of the code). 
Although this is welcome, it should be noted that most implemen- 
tations contain other 'artificial' parts (such as slope limiters used in 
the interpolations), and any user of a code must bear these in mind. 

4. 1.3.1 The AMRA code AMRA is an AMR code developed 
bv lPlewa and Miilleil 1200 ll) . For the disk-planet interaction prob- 
lem we used the Herakles solver which is an implementation of 
the PPM algorithm . HERAKLES was derived from PROMETHEUS 
jFrvxell et alll989h and provides all the functionality of its prede- 
cessor. The original Riemann solver for complex equations of state 
was replaced by a much simple r non-iterative (but still exact) ver- 
sion suited for isothermal flows (Balsara 1994). All problems were 
computed with Courant number of 0.8. The planet was placed in the 
corner of a grid cell, to make the grid layout around it as symmetric 
as possible. 

4.1.3. 2 The FLASH code The FLASH code JFrvxell et alJ 
2000) is an AMR code implementing the PPM algorithm in its 
Direct Eulerian form. 6 The Riemann solver was ported from the 
AMRA code. Two sets of results used FLASH, and we shall refer 
to these as FLASH-AG and FLASH-AP 

The FLASH-AG code was based on release 2.3 of FLASH. 
We patched the code to work as accurately as possible in polar co- 
ordinates, particularly enforcing the conservative transport of an- 

6 FLASH is available at http://www.flash.uchicago.edu/ 



gular momentum. The Courant number was 0.8 in the simulations 
presented here. 

Instead of running in polar co-ordinates, the FLASH-AP 
version of the code used the original Cartesian formulation of 
FLASH. The grid cells were sized to give the same radial res- 
olution, although since the grid went to r = 0, the grid size had 
to be larger than in the cylindrical schemes to achieve the same 
resolution. The code was run at resolution n x x n y = (320, 320) 
and n x xn y = (640, 640). The boundaries were open and there was 
free gas flow inside 0.4a. The damping condition described in Sec- 
tion !3.2l was applied on the outer boundary ring but not in the inner 
disc. The Cartesian grid was fixed in space, and the planet and star 
were free to move over it (integrated using a simple Runge-Kutta 
method). A Courant number of 0.7 was used in the simulations. 

4.1.3.3 The Rodeo code T his code us es the approximate 
Riemann solver suggested b y iRod (l981), and extended by 
Eulderink and Mellema ( 1995) to general non-inertial, curvi-linear 
coordinate systems. The limiter function used is "superbee", and 
unlike PPM-type approaches, limits the characteristic variables, 
rather than the primitive variables. The code uses an AMR scheme 
similar to PARAMESH (used in the FLASH code). The source 
terms are hand led through the so-call ed stationary extrapolation 
method (Euld erink and Mellemalll995l) . which ensures that phys- 
ically stationary solutions remain stationary. The equation of 
state was strictly isothermal. A f ull description can be found in 
iPaardekooper and MellemaH200(j). 

4.1.3.4 The JUPITER code The JUPITER code is a nested 
grid Godunov code, that can be used in Cartesian, cylindrical or 
spherical geometry, either in ID, 2D or 3D. The JUPITER code 
uses a 'two shock' Riemann solver, which a ssumes th at the two 
waves leaving the interface are Shockwaves (Toro 1999) (the code 
can also use a 'two rarefaction' solver, or a full iterative one). The 
rest of the Riemann solver (the sampling of the Riemann fan) is ex- 
act. Assuming that the two waves are Shockwaves is not so bad as it 
might first appear. Firstly, some initial Riemann states do give rise 
to two Shockwaves. Secondly, the differences from the full Rie- 
mann solution are relatively small, so long as the contrast across 
the interface is not extreme. In extra tests (not included here), the 
differences between a 'two shock,' 'two rarefaction' and full Rie- 
mann solver were found to be slight for our comparison problem. 
The predictor step (which provides the left and right states of the 
Riemann problem at the zone interface) is a linear piecewise char- 
acteristic method using the monotonized centred slope limite r, and 
which uses a slope splitting technique (Pem ber et alJ ll995). The 
full viscous stress tensor is conservatively implemented in the three 
geometries. No artificial viscosity was required, and the Courant 
number was 0.7. 

4.1.3.5 The TRAMP -PPM code TRAMP-PPM is a La- 
grangian remap 7 PPM (Woodward and Colellalll984h code. It is 
based on the routines prov ided in the VH- 1 p ackage, modified for 
accretion disk simulations (Blondin and Lufkin 1993). The modifi- 
cations involve adding the conservation of angular momentum and 
equations to treat the evolution of internal energy. Here always the 
full Riemann problem is solved iteratively and we approximate the 



7 Cell boundaries are allowed to move during the advection step, and the 
results are then interpolated back onto the fixed grid 
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isothermal case with y — 1 + 10 -10 . In the current isothermal sim- 
ulations the PPM code does not use any artificial viscosity. This 
implementation works in the corotating frame where the star is the 
centre of the coordinate system, and hence explicitly incorporates 
the extra acceleration terms due to the offset from the centre of 
mass. 

4.2 Particle Based Codes 

Rather than trying to solve the equations of hydrodynamics on 
a grid, a second group of codes decompose a fluid into small 
packets of mass (particles), and then follow their evolution. The 
method in most common use today is that of Smo othed Particle 
Hydr odynamics (SPH), deve l oped independently by iLucvl Jl977l) 
and iGingold and Monaghari ll977r) . We shall describe the basic 
characteristics of SP H now . For a more detai l ed tre atment, the 
reader should consult iBenzl (1990); Monaghan (1992) and refer- 
ences therein and thereto. 

In SPH, each particle's properties are spread (or smoothed) 
over a small volume of space contained within a smoothing length, 
h. For example, smoothing out the particle's mass gives its contri- 
bution to the density at each point in space. The smoothing function 
(or kernel), W(r, h), is not constant, but increases towards the parti- 
cle's position (assumed to be r = 0). In the limit h —> 0, W becomes 
a c) -function, and perfect fluid behaviour is obtained (with an infi- 
nite number of particles). A Gaussian would be a possible choice 
for W, but compact kernels (where W = for r greater than some 
''max) are preferred for computational simplicity. Particles within 
the range of the compact kernel are called the neighbours. When 
appropriate to the problem, modern SPH codes will allow each par- 
ticle to have its own smoothing length, chosen to keep the number 
of neighbours constant (typically a few tens). The smoothing length 
is also used to limit the timestep in a way similar to the CFL con- 
dition mentioned above. 

The major advantage of SPH is that its particle nature makes 
it fully Lagrangian: there are no advective terms in the equations 
of motion. This makes the codes more straightforward to write 
and understand. Since high densities imply that more particles are 
present, 8 SPH naturally concentrates resolution in high density re- 
gions. Good use can be ma de of this in collapse simulations (e.g. 
iDelgado-Donate et al|2003|). 

However, there are disadvantages too. Foremost is the matter 
of viscosity. SPH requires an artificial viscosity to prevent inter- 
particle penetration, and this tends to make SPH codes quite dis- 
sipative. Resolution can also be a problem in certain calculations. 
For example, in the disc calculations presented here, most of the 
particles are going to be in the outer portions of the disc, and not 
doing very much. Also, the details of the gap are of most interest, 
and SPH will have fewer particles there. 

4.2.1 The SPHTREE code 

This code owes its name to the tree used to locate particle neigh- 
bours. The calculations presented here used 250 000 particles for 
the disc, with the star and planet being point masses. SPH particles 
that move to within an accret ion radius of either the star or planet 
are accreted feate et all 19 95) but in the case of the planet, once the 
initial ramp up is complete, we do not allow its mass to increase. We 

8 Although it is possible to let particle masses vary in SPH, it is not entirely 
trivial to do so 



use the standard SPH viscosity (e.g. Monaghar i 1992h. with a = 0.1 
and (5 = 0.2, but also can include the Balsara switch fealsarall 995) 
to reduce the shear component of the artificial viscosity (see also 
lLodato and Ricel2004h . A huge saving in computat ional time is ob- 
tained by using individual particle timesteps ( Bate et al. 1995) with 
the time-steps for each particle lim ited by the Courant condition, 
a force condition (Monaghan 1992) and a Runge-Kutta integrator 
accuracy condition. 

4.2.2 The PARASPH code 

ParaSPH is a parallelized (using MPI) smooth particle hy drody- 
namics code. It follows the approach of Fleb be et alJ dl994l) . solv- 
ing the Navier-Stokes equation including the entire viscous stress 
tensor. In contrast to the usual approach of an artificial viscosity of 
Monaghan and Gingold ( 1983), we use an artificial bulk viscosity. 
This allows for an accurate treatment of the physical shear viscos- 
ity and for easy comparison to the grid code results, since a con- 
stant kinematic viscosity coefficient can be modeled. Additionally 
we use the XSPH device to pr event particles from mutual penetra- 
tion (see e.g. lMonagha n 1989) Variable smoothing lengths keep the 
number of neighbours at 75. The time integration is performed us- 
ing a fourth order Runge-Kutta-Cash-Karp integrator for both the 
particles and the pl anet. The code is described in more detail in 
Schaf er etallbOQA . 

We do not implement exactly the boundary conditions de- 
scribed in sect. 13.21 Instead we add virtual particles to the simu- 
lation. They are assigned all physical relevant quantities, such as 
density, velocity and so on, but are kept in Keplerian orbit about the 
star. By their interactions, the virtual particles prevent the SPH par- 
ticles from escaping. For the calculations presented here, we used 
300 000 SPH particles and 50 000 boundary particles. 



5 RESULTS 

In this section we present the results for each of the runs. The sim- 
ulations are run for up to 500 orbital periods using the codes de- 
scribed in Section |4] We compare the contours of surface density, 
vortensity and averaged density profiles obtained in the numerical 
calculations at several times. The time evolution of the grid mass 
and gravitational torque acting on the planet are shown divided in 
several contributions. The Fourier transform of the torques is cal- 
culated to investigate the influence of vortices and disc eccentricity 
on the torque acting on the planet. Several basic properties of the 
disc-planet system are discussed based on the agreement between 
the codes. In Section 1531 we study how the difference between the 
codes change as the numerical resolution increases. 

The comparative surface density and vortensity maps are 
shown for each scheme in the order they appear in Section |4] Note 
that TRAMP-PPM and TRAMP- vanLeer were only run for the 
inviscid setups, while SPHTREE and PENCIL were run for the vis- 
cous cases (see Table 0. Figure [I] shows the legend used in the 
surface density profiles, mass and torque evolution plots in this 
Section. Different types of algorithms are plotted with different 
linestyles. 

5.1 Inviscid Jupiter 

Firstly we consider the case of a Jupiter embedded in an inviscid 
disc. The planet fixed at a given radi us opens a deep gap in the 
disc as predicted by standard theory iLin and Papaloizou 1986a; 




Figure 2. Density contours in logarithmic scale after 100 orbits for the inviscid Jupiter simulations with overplotted theoretical prediction of the planetary 
wake position. The codes are presented in the same order as in Section [4] The upwind methods' results are displayed in the first panels followed by the 
shock-capturing codes and lastly the particle based codes. The density scale ranges between —1.7 < log(£/Eo) < 1- 
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NIRVANA-GDA 

- NIRVANA-GD 
NIRVANA-PC 
RH2D 

- GLOBAL 
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■ Pencil 

AMRA 
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Flash-AP 
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- JUPITER 

TRAMP-PPM 

SPHTREE 

ParaSPH 

Figure 1. Common legend for the comparative plots in section [5] Upwind 
codes are represented by solid lines, shock-capturing codes by dotted lines 
and SPH codes by dashed lines. 

Ward and Hahn l200(t) . The contrast between the initial density and 
the deepest regions in the gap is about 2 orders of magnitude af- 
ter 100 orbits. The planet forms strong trailing spiral arms due 
to the differential Keplerian rotation close to the Lindblad reso- 
nances. Low density regions start to develop behind the shocks 
where the fluid elements encounter the shock at high pitch an- 
gle and change their trajectory. These regions travel in horseshoe 
orbits in the corotation region cl earing the gap as described in 
Korvcanskv and PaDaloizou ( 1996) and creating locked fluid areas 
at the Lagrangian points inside the gap. 

Figure 13 shows the density contours at 100 orbits for all the 
codes. The dashed line represents the theoretical p osition of the 
shock wave predicted bv lQgilvie and Lub ow (2002). In the Jupiter 
simulations, the planet mass is too high for this theoretical estima- 
tion but it allows us to compare the spiral arms pitch angle between 
the models. The planetary wakes have a high pitch angle compared 
with the theoretical calculation next to the planet. There is a sec- 
ondary shock in the Eulerian codes which starts near the L5 point 
and has approximately the same opening angle as the theoretical 
prediction. The secondary shock seems to be related with the den- 
sity excess inside the gap behind the planet. In the outer part of the 
disc, the pitch angle of the primary and secondary shocks are very 
similar. The existence of secondary shocks and the tight ness of th e 
spirals depend primarily on the equation of state used lKlevll 999). 

There are two density enhancements in all grid-based models 
located close to the L4 and L5 points at azimuthal distance A0 = 
±7r/3 from the planet. In the SPH and FLASH-AP codes the gap 
is almost completely clean. Theoretically, the calculation should 
produce a nearly symmetric density distribution inside the gap at 
both sides of the planet's location for the case of a planet in a fixed 
orbit, which is observed in our results. 

Shock-capturing codes that use cylindrical coordinates such as 
FLASH-AG, AMRA and TRAMP-PPM have filamentary struc- 
ture visible in the disc and the gap possibly due to the high-order 
scheme of the codes. The filaments can be produced by instabil- 



ities generated locally on the corrugated spiral shock. Their an- 
gle does not match the angle of the spiral shocks, so they can- 
not be generated around the planet. In tests performed with the 
FLASH -AG code, the filaments appear in high resolution calcu- 
lations with larger amplitude but the same structure. The shock- 
capturing codes using a different algorithm than PPM such as 
JUPITER and Rodeo do not present filaments although they seem 
to have more structure in the disc than the upwind methods' results. 

Figure|3] shows the vortensity contours calculated in the coro- 
tating frame for the different models. There are bumps rotating 
along the edges of the gap opened by the planet in the grid codes 
in cylindrical coordinates which survive until the end of the sim- 
ulations. The resolution does not permit us to determine whether 
these density lumps have locally rotating flow around the core of 
the vortex. The vortices are larger in the upwind schemes. Af- 
ter 100 periods, most codes show a single bump rotating along 
the outer edge, although NlRVANA-GD, AMRA, FLASH-AG 
and TRAMP-PPM have two bumps which eventually merge by 
200 periods. The knots are dominant in the AMRA, RODEO and 
TRAMP-PPM simulations and generate their own spiral shocks 
which extend into the disc. Most of the codes show one or several 
smaller density excesses at the inner edge. The vortices in the outer 
disc interact with the planetary shock and generate quasi-periodic 
oscillations in the spiral arms. The oscillations could also be pro- 
duced by instabilities near the planet that interact with the blobs 
moving along the edge of the gap and are propagated along the 
shocks. Reflected waves appear in the NlRVANA-GDA, RH2D, 
GLOBAL and GENESIS codes (see Figure |3|, despite the use 
of wave killing boundaries. 

In the PARASPH code the gap edges are less steep than in 
the case of the grid based calculations possibly due to the artificial 
viscosity. The planetary wake is weak and almost not visible in the 
inner disc. 

The azimuthally averaged density profiles and their residuals 
normalized by the disc mass after 100 orbits are plotted in Figure^! 
Note that the relevant portion of the domain is that between 0.5 and 
2.1 a, since wave killing conditions are implemented next to the in- 
ner and outer boundaries. The depth and width of the gap is in good 
agreement for the grid based models, with a slightly wider gap for 
the AMRA code. FLASH-AP has a more depleted inner disc due 
to the open inner boundary condition implemented in Cartesian co- 
ordinates. The Cartesian geometry is also more diffusive in this 
problem and has a lower resolution close to the primary compared 
with the polar coordinates codes. On the other hand, the shape of 
the outer disc in FLASH-AP's profiles is similar to the cylindri- 
cal codes profiles. A wider gap is seen in the PARASPH simula- 
tion. The oscillations seen in the outer disc are also consistent in 
all the codes with a local maximum in the FLASH-AP profile at 
2a. The density peaks close to the edges of the gap — specially in 
the inner disk — have a larger spread which is associated with the 
size of the vortices. Shock-capturing codes have smaller vortices 
in the outer edge than the upwind methods. The maximum at the 
planet location is higher on average for the shock-capturing codes. 
The PARASPH code has smoother profiles farther away from the 
planet position due to the fact that the planetary wakes are smeared 
out. The residuals of the averaged profiles divided by the disc mass 
with respect to the mean value are shown in the bottom panel in 
Figure |4] Since the total disc mass is different after 100 orbits for 
the various models, the density profiles normalized by the disc mass 
have in general a better agreement. However, the PARASPH and 
FLASH-AP codes have most of the mass loss in the inner disc 




Figure 3. Vortensity contours in logarithmic scale after 100 orbits for the inviscid Jupiter models. The vortensity range is —0.5 < log(Q < 2. 
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Figure 4. The upper panel shows the normalized surface density profiles 
averaged azimuthally over 2k after 100 orbits for the inviscid Jupiter runs. 
In the lower panel, the differences between each model and the mean value 
are shown in logarithmic scale as log(E/Mdi sc ) - (log(E/Mdi sc )), where the 
angle brackets represent the mean. The surface density has been divided by 
the disc mass at 100 periods to remove the dependence on the mass loss due 
to the boundary conditions. 



and this method may artificially increase their residuals in the outer 
disc. 

We plot the density slices opposite to the planet after 100 or- 
bits in Figure [5] The width and depth of the gap agree well for the 
different codes but with a larger dispersion than in the averaged 
profile. The amplitude of the peaks at the edges of the gap differ 
since there are vortices which have different sizes and positions 
with respect to the planet at a given time. On the other hand, the 
size and the position of the wave crests agree within a few percent. 

In Figure [6| the azimuthal cuts of the surface density maps at 
the planet radius are displayed. There is a sharp density spike at 
the planet position. The shape and prosition of the density bumps 
at L4 and L5 is slightly asymmetric. The peak at the trailing La- 
grangian point L5 is larger than at the leading L4 point for all Eu- 
lerian codes, with more conspicuous peaks and larger asymmetry 
in the shock-capturing schemes. In the FLASH-AP results there 
are asymmetric bumps at the Lagrangian points in the beginning 
of the simulation, but they have disappeared at 100 periods. In the 
PARAS PH calculation the gap is almost completely cleared and 
no bumps at the Lagrangian points are observed. PARAS PH has 
also a smaller peak at the planet location. 

The disc mass loss rate evolution is plotted in Figure for 
the Eulerian codes. The total disc mass is not conserved due to the 
wave damping condition described in section IT2I There is a larger 
mass loss rate in the FLASH-AP code owing to the mass accre- 
tion in the inner disc but it reaches an equilibrium value consistent 
with the cylindrical codes at late times. Some codes gain mass at 




Figure 5. Surface density profiles opposite to the planet position after 100 
orbits for the inviscid Jupiter runs. 
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Figure 6. Surface density azimuthal slice at the planet radius after 100 or- 
bits for the inviscid Jupiter calculations. The trailing Lagrangian point L5 is 
located at azimuth ~ — 1/3 and the leading Lagrangian point L4 is at ~ 1/3 
in the normalized azimuthal units. 



the beginning of the simulation and start losing mass after about 
10 orbits. The total mass after 200 orbits is reduced by about 8% 
in the AMRA and FLASH-AG codes which use shock-capturing 
algorithms. Other codes like NIRVANA- GD A, NIRVANA- GD and 
RH2D show a smaller mass loss of about 3%. The outer disc mass 
decreases slowly and in some codes like JUPITER it remains al- 
most constant during 200 orbits. During the first few orbits, when 
the gap is not completely cleared, there is material flowing from 
the outer to the inner disc perhaps due to the artificial viscosity. 
The inner disc mass shows a strong decrease, especially in the 
shock-capturing codes. Despite the spread in mass loss for different 
codes, the surface density do not show strong variations between 
the codes. 

The waves excited by the planet deposit angular momentum in 
the disc when they are dissipated. There is an initial smooth phase 
where the torque increases in absolute value during the first few 
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Table 5. Averaged gravitational torques between 175-200 periods in units 
where a = 1, P = 2n and M* = 1 — /I for the Jupiter inviscid simulations at 
the end of the simulations. The time is given in orbital periods of the planet. 



Code Torque 



Nirvana- GD A 


-1.452354 x 10" 


-5 


Nirvana-PC 


-1.512417 x 10" 


-5 


RH2D 


-1.930871 x 10" 


-5 


GLOBAL 


-1.550768 x 10" 


-5 


GENESIS 


-1.565123 x 10" 


-5 


TRAMP-vanLeer 


-1.818716 x 10" 


-5 


AMRA 


-3.769203 x 10" 


-5 


FLASH-AG 


-7.014221 x 10" 


-5 


FLASH-AP 


-2.187462 x 10" 


-5 


Rodeo 


-1.880762 x 10" 


-5 


JUPITER 


-2.721250 x 10" 


-5 


TRAMP-PPM 


5.611237 x 10" 


6 



Figure 7. Evolution of the disc mass loss rate over 200 orbital periods for 
the inviscid Jupiter simulations. 
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Figure 8. Profiles of total specific gravitational torque as a function of ra- 
dius after 100 orbits for the inviscid Jupiter simulations. 



orbits while the planet is growing. Afterwards, the torques start to 
display strong oscillations at the time when the vortices are cre- 
ated. Vortices grow due to the steep gradients at the gap edges 
and through interaction with the planetary wakes. We do not have 
enough time resolution in the density snapshots to follow the vor- 
tex formation and evolution. The mean value decreases slightly in 
most codes until the point when the gap is completely cleared and 
stays roughly constant up to the end of the simulation. The effect 
of the large torque oscillations on the planet migration needs to be 
studied with a free moving planet. 

The torque from within the Roche lobe show significant dif- 
ferences between the codes. The density has a local maximum at 
the planet location which depends on the interpolation order of the 
code, although the total mass inside the Roche lobe is similar. The 
planet is not located in a cell's corner in all codes and this causes 
asymmetries in the mass distribution around the planet In addition, 
the region close to the planet is not well resolved at our resolution. 



In the following discussion, we compare the torques excluding the 
contribution from the Roche lobe. 

Figure|8]shows the profiles of the derivative of the total torque 
excluding the Hill sphere with respect to the radius. The time de- 
pendence of the vortices position with respect to the planet pro- 
duces a rapidly changing torque. Therefore, the different codes 
have different specific torque profiles at a given time. The varia- 
tion appears close to the gap edges where most of the angular mo- 
mentum is deposited. The differences are larger at the outer edge 
position where the vortices are bigger than at the inner edge. On 
the other hand, farther away from the planet position the torques 
are remarkably similar for all the codes. 

The time evolution of the gravitational torque acting on the 
planet is shown in Figure [9| divided in inner, outer disc and to- 
tal contributions. A running time average over 10 orbital periods 
has been applied to the data to avoid large oscillations. The torque 
contribution from the disc material inside the planet's orbit gives a 
positive torque on the planet which tends to drive the planet out- 
wards in all models, while the torque from the material outside 
the planet's orbit pushes the planet towards the star. The outer disc 
contribution is dominant and gives a total negative torque which 
takes away angular momentum from the planet and would cause 
inwards migration in case the planet were released. The torque or- 
der of magnitude and sign agrees for all the codes after 200 orbits 
except for TRAMP-PPM which has a value close to zero. The av- 
eraged values at the end of the simulation are given in Table|5] 

The power density spectra (PDSs) of the corresponding grav- 
itational torque components are shown on the right hand side plots 
in Figure |9] The panels show the low frequency part of the power 
spectrum in logarithmic scale. The semi-periodic oscillations asso- 
ciated with vortices rotating along the gap edges are present in the 
PDSs for models where blobs appear in the density maps next to 
the gap. Several peaks appear in the outer disc PDS with frequency 
corresponding to roughly 0.4 times the planet's orbital frequency 
with several harmonics. This frequency is the difference between 
the planet's orbital frequency and the angular velocity of the high 
density regions moving next to the gap. Assuming that the density 
lumps orbit the central star with Keplerian speed, the position of 
the blob estimated from the PDS frequency is about 1.4a, in agree- 
ment with the center of the blobs observed in the density maps. The 
harmonics possibly appear because the potential of an extended 
density blob is not sinusoidal and creates amplified multiple fre- 
quencies. In some codes, there are several vortices next to the outer 
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Figure 9. Torques time series for the inviscid Jupiter simulations smoothed over 10 periods and the corresponding normalized power density spectra in 
logarithmic scale. The upper panels show the torque contribution from the inner disc, the middle panels show the torque from the outer disc and in the bottom 
panels the total torque is plotted. In all plots the contribution from inside the Hill sphere is ignored to avoid numerical noise. 



edge which perturb the planet with the same frequency but different 
phase. In the inner disc contribution PDS, there are several codes 
with peaks at about 0.7 times the planet's orbital frequency and its 
harmonics. The estimated blob position is about 0.7a, which again 



agrees with the center of the vortices observed next to the inner 
edge of the gap in the density maps. 

The PDSs of the torque from the material inside the Roche 
lobe show high frequency quasi-periodic variations at about 20 
times the planet orbital frequency. This high frequency oscillations 
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Figure 12. The upper panel shows the surface density profiles averaged 
azimuthally over the whole azimuthal range after 100 orbits for the viscous 
Jupiter case. In the lower panel, the difference between each model and the 
mean value is shown as defined in Figure^ 



may be caused by the circumplanetary disc which makes several or- 
bits around the planet within the planet orbital period although the 
region is poorly resolved at our resolution. There is a local max- 
imum in density inside the Roche lobe and the material gives a 
leading contribution to the total torque acting on the planet. 



5.2 Viscous Jupiter 

The density contours for Jupiter in a disc with Navier-Stokes vis- 
cosity v = 10 -5 are shown after 100 orbits for all the codes in Fig- 
urelTUl The planet opens a narrower gap in the disc in this case. The 
flow is much smoother than in the inviscid calculation and the blobs 
moving along the gap are not observed. The density enhancements 
seen at the Lagrangian points inside the gap in the inviscid calcu- 
lations are not present. The spiral arms generated by the planet are 
stationary. The filamentary structure that appeared in the inviscid 
Jupiter runs in the shock-capturing codes is reduced in amplitude. 
The reduction is stronger in FLASH -AG and FLASH-AP results 
than in AMRA which uses a different dissipation algorithm. 

In Figure [H] we plot the vortensity for the viscous Jupiter 
case. The maps are smooth compared with the inviscid simula- 
tions and vortices are not visible in the logarithmic scale. Reflected 
waves are visible in the NlRVANA-GDA, RH2D, GLOBAL and 
GENESIS results despite the use of the wave killing zones. 

In Figure [l2l we show the azimuthally averaged density pro- 
files and normalized residuals after 100 orbits. The depth and width 
of the gap agree well for the grid codes with a shallower gap in the 
FLASH-AP code. The gap is wider and deeper in the PARASPH 
simulation. The SPHTREE code has a small peak at the planet ra- 
dius and the inner disc is depleted due to mass loss. An slightly 




Figure 13. Surface density profiles opposite to the planet position after 100 
orbits for the viscous Jupiter runs. 
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Figure 14. Surface density azimuthal cut at the planet position after 100 
orbits for the viscous Jupiter runs. 



asymmetric gap structure is observed in most codes with a deeper 
dent outside the planet's orbit. The oscillations in the outer disc 
have disappeared in the grid codes or have been reduced consider- 
ably by the viscosity. The differences of the averaged profiles with 
respect to the mean value are shown in FigurefT2l 

We plot the surface density profiles at = n after 100 orbits 
in Figure [H] The peaks at the edges of the gap agree well since 
there are no vortices in the viscous runs and the spiral arms are 
weaker. Due to the viscosity, the gap is narrower and shallower 
than in the inviscid case. The shape of the spiral arms agree within 
a few percent for the grid based codes. The SPH codes agree in the 
general shape of the density profile but have weaker spiral waves. 
SPHTREE has a density peak in the middle of the gap opposite 
from the planet. 

In Figure Q4| we plot the azimuthal cuts of the surface density 
maps at the centre of the gap. A sharp density spike is seen at the 
planet position in all codes. The density bumps at the equilibrium 
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Figure 10. Density contours after 100 orbits for the viscous Jupiter simulations. The dashed line is the estimated theoretical position of the planetary shocks. 
The density range is again —1.7 < log(L/Eo) < 1- 
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Figure 11. Vortensity maps for the viscous Jupiter case after 100 orbits. The logarithmic scale is -0.5 < log (£) < 2. 
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Figure 15. Evolution of the disc mass loss rate over 200 orbital periods for 
the viscous Jupiter case. 



points inside the gap have disappeared. Most of the grid codes show 
a constant density of about 15% of the initial value. PARASPH 
has a lower density than the grid codes, while SPHTREE shows 
a density about twice as large as the grid codes. The presence of 
oscillations in the SPHTREE azimuthal profile may be explained 
because the number of particles is too small to resolve the gap. The 
effective resolution after projection in the radial range [0.4a, 2.5a] 
is smaller than in the grid simulations since the radial domain ex- 
tends until 10a and at the end of the simulation a significant fraction 
of the particles has been accreted. 

We plot the evolution of the disc mass loss rate in Figure fT5l 
There is less mass loss than in the inviscid Jupiter case due to the 
weaker waves and the agreement in the loss rate is generally very 
good. FLASH-AP has a larger mass decrease due to the open in- 
ner boundary. Pencil has a very small mass loss possibly due to 
the freezing zones in the boundaries. The total mass loss after 200 
orbits shows better agreement than in the inviscid case. Upwind 
methods show a reduction of about 5% of the initial mass, while 
RODEO, AMRA and FLASH- AG codes lose close to 8% of their 
mass. During the first few orbits, there is again gas flow from the 
outer to the inner disc when the gap is not cleared. The outer disc 
mass decreases slightly for some codes while others present an in- 
crease of roughly 1%. There is a substantial decrease in the inner 
disc mass with an agreement of approximately 10% between the 
different models. 

The amplitude of the torque oscillations are smaller compared 
with the inviscid runs. There is again an initial stage where the 
torque increases in absolute value while the planet mass is increas- 
ing. The torques start to oscillate at about 10 orbits and later possi- 
bly due to the formation of small vortices or eccentricity of the disc. 
In most codes the oscillations decrease and become very small by 
the end of the simulation. 

In Figure ITol the profiles of the specific total torque exclud- 
ing the Hill sphere are shown. The profiles show a much better 
agreement than in the inviscid Jupiter simulations. In the viscous 
case, vortices are not observed in the density maps after 100 orbits 
and the torque radial profiles are not time dependent. There is a 
dominant contribution from the corotating region in the grid-based 
schemes from the exchange of angular momentum with gas flowing 
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Figure 16. Profiles of total specific torque acting on the planet after 100 
orbits for the viscous Jupiter case. 

in horseshoe orbits, although the material inside the Roche lobe is 
not considered. The outer disc gives a negative torque contribution 
on the planet driving inwards migration and the inner disc produces 
a positive torque that pushes the planet outwards. The profiles of the 
polar coordinates codes agree within a few percent. 

The time average of the torque acting on the planet and their 
PDSs are shown in Figure [H] The outer disc torque contribution 
is again dominant and gives a negative total torque. The total aver- 
aged torques at the end of the simulation are shown in Table [6| 
The PDSs of the different torque contribution are shown in the 
right hand side panels in Figure [n] The plots show the low fre- 
quency part of the PDSs in logarithmic scale. There is a peak at 
0.3 times the planet's orbital frequency and several multiples in the 
outer disc PDS. In some models, there is also a small peak at the 
same frequency in the PDS form the inner disc. This quasi-periodic 
oscillations may be produced by vortices appearing during the first 
orbits of the simulation and eventually removed by the viscosity. 
Other possible explanations are asymmetry in the edge of the gap 
or slight eccentricity of the disc. 

The torque from the gas inside the Hill sphere presents again 
a power spectrum with high frequency peaks at several times the 
Keplerian frequency at the planet radius. The smoothing length is 
close to half of the Hill radius and the resolution in the Roche lobe 
is low to study the possible presence of a circumplanetary disc ro- 
tating at high angular frequency. 



5.3 Inviscid Neptune 

The dip opened by Neptune after 100 orbital periods is much shal- 
lower than for the Jupiter case. The surface density maps are plot- 
ted in FigurelT8lfor a Neptune mass planet embedded in an inviscid 
disc. The spiral arms created by the planet are significantly weaker 
than in the Jupiter calculations and are in better agreement with the 
theoretical prediction of the shock positions shown by the dashed 
line. In the SPH simulations the shocks are extremely weak. There 
are no overdense regions around the Lagrangian points inside the 
gap in any of the calculations since the gap is not deep enough. 
Along the edge of the gap there are several blobs in the AMRA, 
RODEO and TRAMP-PPM results, which are smaller than in the 
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Figure 17. Running time averaged torques for the viscous Jupiter simulations and the corresponding PDSs of the raw data. The plots are shown in the same 
order as in Figure |9| All the figures exclude the Roche lobe contribution. 



inviscid Jupiter calculations. The FLASH, AMRA and JUPITER 
codes show ripples in the disc and the gap with lower amplitude 
than in the inviscid Jupiter simulations. 

The comparative vortensity maps in the corotating frame are 
shown in Figure [l9| Vortices moving along the gap are observed in 
the grid codes, although they are smaller than in the Jupiter case. 



The azimuthally averaged density profiles after 100 orbits are 
plotted in Figure |20| The depth and width of the gap is again in 
fairly good agreement for the Eulerian codes. FLASH- AG's gap is 
shallower than the other grid based codes. PARAS PH has a wider 
and deeper gap than the grid models and a depleted inner disc. The 
gap profile of the Eulerian codes is slightly asymmetric with the 
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Figure 18. Surface density contours after 100 orbits for the inviscid Neptune simulations. The theoretical estimation of the spiral wakes is represented by the 
dashed line. The density scale ranges between -0.4 < log(L/Eo) < 0-3- 
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Figure 19. Vortensity contours in logarithmic scale after 100 orbits for the inviscid Neptune calculations. The vortensity range is —0.1 < log(£) < 1. 
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Table 6. Averaged torques between 175-200 periods in units where a = 1, 
P = 2k and M* = 1 — ji for the Jupiter viscous simulations. 



Code Torque 



Nirvana- GD A 


-7.279010 x 


10 


-5 


Nirvana-PC 


-8.591188 x 


10 


-5 


RH2D 


-7.895477 x 


10 


-5 


GLOBAL 


-5.424962 x 


10 


-5 


GENESIS 


-7.980679 x 


10 


-5 


Pencil 


- 1.265820 x 


10 


-4 


AMRA 


-7.269365 x 


10 


-5 


FLASH-AG 


-9.288739 x 


10 


-5 


FLASH-AP 


-6.681127 x 


10 


-5 


Rodeo 


-1.061018 x 


10 


-4 


JUPITER 


-8.067247 x 


10 


-5 
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Figure 20. Surface density profiles averaged azimuthally over In after 100 
orbits for the inviscid Neptune runs are shown in the upper panel. The resid- 
uals in the lower panel are defined as in Figure|4| 



deepest part just outside of the planet radius. In the lower panel in 
Figure |20| we show the residuals of the averaged profiles divided 
by the disc mass with respect to the mean value. 

In Figure E] we plot the surface density at = %. The shape 
and amplitude of the waves in the disc agree well for the different 
codes outside the wave damping boundaries. There is a larger dis- 
persion at the inner gap edge and in the middle of the gap for the 
shock-capturing schemes. The gap is slightly asymmetric for the 
majority of the codes. 

In Figure l22l the azimuthal slices of surface density at the 
planet position are shown. A large density peak is observed again at 
the planet position for all the grid codes. The FLASH- AG, RODEO 
and JUPITER density in the centre of the gap after 100 orbits 
is close to the initial density with depressions next to the planet. 



Figure 21. Surface density profiles opposite to the planet position after 100 
orbits for the inviscid Neptune runs. 
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Figure 22. Surface density azimuthal slice at the planet location after 100 
orbits for the inviscid Neptune calculations. 



The rest of grid codes show a density decrease of about 10-20%. 
PARAS PH has the lowest density inside the gap with a decrease 
of about 60% from the initial value. 

Figure [23] shows the grid mass loss rate as a function of time 
for the grid-based codes. All models show total mass loss due to the 
wave killing condition. The FLASH-AP code has mass loss in the 
inner disc due to the absence of a solid inner boundary in Cartesian 
coordinates but converges to a value of a few times 10 -5 after 200 
orbits. There is mass increase in the inner disc for some schemes in 
the beginning of the simulation. This suggests that there is gas flow 
trough the gap from the outer to the inner disc since in the Neptune 
simulations the gap is shallower and the planet generates weaker 
shocks. The artificial viscosity may cause the flow from the outer 
to the inner disc. Another possible explanation is that the damping 
wave condition near the inner boundary adds artificially angular 
momentum to the disc. 

We plot the profiles of the derivative of the torque with re- 
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Table 7. Averaged torques at the end of the simulations in units where a - 
P — 2k and M* — 1 — /I for the Neptune inviscid simulations. 
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Figure 23. Disc mass loss rate evolution for the inviscid Neptune simula- 
tions. 
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Figure 24. Specific torque profiles after 100 orbits for the inviscid Neptune 
simulations. 



spect to the radius in Figure |24| The agreement between the codes 
is good compared with the inviscid Jupiter case, especially for the 
cylindrical grid hydrodynamical codes. The vortices do not perturb 
the planet strongly and the specific torque radial profiles are sta- 
tionary. The outer disc generates again a negative torque acting on 
the planet and the inner disc gives a positive torque. At several Hill 
radii away from the planet location the torques become negligible. 

In Figure[25] the time average of the gravitational torques act- 
ing on the planet and their associated PDSs are plotted. The total 
torque after 200 periods agree within a factor 2 (see Table [7J. up- 
wind results show good agreement while the shock-capturing re- 
sults have larger oscillations. The oscillations observed in the raw 
data and PDSs may be produced by short-lived vortices appearing 
during the first orbits which are not visible at later time in the den- 
sity snapshots. 




2.5 



Figure 28. The upper panel shows the surface density profiles averaged az- 
imuthally over 2k after 100 orbits for the viscous Neptune runs. The resid- 
uals in the lower panel are defined as in Figure|4| 



5.4 Viscous Neptune 

In Figure|26]the comparative surface density contours after 100 or- 
bits for the viscous Neptune case are pl otted. The theoretical esti- 
mation of the spiral shocks positions bv lQgilvie and LubowH 2002) 
is shown by the dashed line. The flow is smoother than in the in- 
viscid Neptune simulations. The density lumps moving along the 
edge of the gap have disappeared and the planetary wakes are sta- 
ble. The filamentary structures in the shock-capturing simulations 
have a reduced amplitude compared with the inviscid case. 

The vortensity maps are shown in Figure l27l The density blobs 
lying next to the gap's edge are not observed in the viscous simu- 
lations in logarithmic scale. Several codes show wave reflection at 
the outer boundary despite the wave damping condition. 
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Figure 25. Time averaged torques and their PDSs for the inviscid Neptune simulations. The plots are shown in the same order as in Figure|9|and exclude the 
material inside the Roche lobe. 



The smoothed density profiles are shown in Figure l28l for 
the viscous Neptune calculations. The gap profile is again in good 
agreement for the polar grid hydrodynamics codes. The gap is shal- 
lower for FLASH- AG than for the other Eulerian codes. FLASH- 
AP has a wider and deeper gap with a flat shape. PARAS PH has 
a very deep gap and SPHTREE has a strong asymmetry with the 
deeper depression outside the planet postion. The residuals of the 



averaged profiles divided by the disc mass are shown in the bottom 
panel in Figure |4| 

The surface density opposite to the planet after 100 orbits 
is shown in Figure l29l FLASH- AG has a shallow gap whereas 
FLASH-AP has a deeper and broader gap. The waves observed in 
the inner and outer disc agree within a few percent for the Eulerian 
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Figure 26. Surface density maps after 100 orbits for the viscous Neptune simulations. The dashed line is the estimated theoretical position of the spiral arms. 
The density range is -0.4 < log(£/Eo) < 0-3- 
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Figure 27. Vortensity contours after 100 orbits for the viscous Neptune simulations. The vortensity range is —0.1 < log(£) < 1. 
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Figure 29. Surface density profiles opposite to the planet position after 100 
orbits for the viscous Neptune case. 
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Figure 30. Surface density azimuthal slice at the planet location after 100 
orbits for viscous Neptune calculations. 

schemes. PARAS PH has a very open gap and SPHTREE has a 
noisy profile with a deep cavity outside the planet's radius. 

The surface density azimuthal slices after 100 orbits are plot- 
ted in Figure |30| A density spike appears again at the planet lo- 
cation in the grid codes. The grid codes show a density decrease 
of approximately 10-20% of the initial density in the center of the 
gap, while PARAS PH has a decrease of about 60% as in the invis- 
cid Neptune case. 

The disc mass loss rate is shown in Figure I3T1 All the models 
apart from FLASH- AP show total mass loss after 200 periods with 
final values consistent within a factor of about 3. Rodeo has a 
sharp jump in mass loss rate at about 155 periods. FLASH-AP 
results have a considerable mass transfer from the region close to 
the star due to the gravitational softening. 

The dT/dr profiles after 100 orbits are shown in Figure l32l 
The profiles show a good agreement between the grid-based codes. 

We plot the time averaged torques acting on the planet on the 
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Figure 31. Disc mass loss rate evolution for the viscous Neptune simula- 
tions. 
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Figure 32. Specific torque profiles after 100 orbits for the viscous Neptune 
simulations. 

left hand side of Figure l33l The torque contribution from the inner 
disc is positive, while the outer disc contribution is negative. The 
outer disc dominates the total torque and would cause an inwards 
orbital shift for a free moving planet. In Table [8] we show the av- 
eraged torques at 200 orbital periods. The torques PDSs are shown 
on the right hand side of Figure [33] The spectrum is rather flat for 
all codes which agrees with the absence of vortices or eccentricity 
in the disc. 



5.5 High resolution simulations 

We studied the convergence of the results running the test prob- 
lem at 2 and 4 times the original linear resolution with some 
of the codes. NlRVANA-PC and NlRVANA-GD Jupiter simula- 
tions were run at resolution n x xn§ = (256,768). Several tests at 
n Y xri0 = (512, 1536) were done with RH2D, NlRVANA-GD and 
FARGO codes for Jupiter and Neptune planet masses. PARAS PH 
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Figure 33. Time averaged torques and corresponding PDSs for the viscous Neptune case. The plots are shown in the same order as in Figure|9]and exclude 
the material inside the Roche lobe. 



was run using 853280 particles and 146720 boundary particles for 
the Jupiter viscous case. 

In the grid-based schemes, the flow is observed to be smoother 
and more stable in time than in the low resolution runs. Vortices are 
still visible in the Jupiter inviscid simulations in the NIRVANA- PC 
and FARGO simulations. The vortices are more extended than in 
the lower resolution calculations and interact with the primary and 



secondary shocks. There is more mass piling up inside the Roche 
lobe after 200 orbit s in the higher resolu tion cases in agreement 
with the results of ID' Angelo et alJ fc005l) . Nevertheless, the aver- 
aged density profiles are very similar to the results presented in the 
previous sections. The gravitational torques are similar in the grid- 
based codes and in good agreement with the low resolution results. 
The PARASPH results with ~ 850000 particles have stronger 
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Table 8. Averaged torques in the window 175-200 periods in units where 
a = 1, P = 2k and M* = 1 — ji for the Neptune viscous simulations. 
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shocks and the density profiles are in good agreement with the grid- 
based results. This suggests that SPH schemes need higher resolu- 
tion to model accurately the corotation region and planetary wakes. 

6 DISCUSSION 

In this paper, we have studied a planet in a fixed orbit embedded in a 
disc using 17 different SPH and Eulerian methods. The codes used 
in the comparison have been thoroughly tested in problems with 
known analytical solutions. The goal of this project was to investi- 
gate the reliability of current astrophysical hydrodynamic codes in 
the disc-planet problem, and to provide a reference for future cal- 
culations. Performing this comparison also aided in the debugging 
of the codes. 

The results show good agreement on the general picture, al- 
though there are some differences in the details. The density maps 
and averaged profiles are consistent for the grid-based methods. 
The variation in the disc mass is of the order of 10% after 100 or- 
bital periods, but this does not seem to produce big differences in 
the surface density distributions. The different boundary conditions 
tested in FARGO do not affect the results since the goal in both 
boundary implementations was to avoid the reflection of waves. A 
preliminary study of convergence using finer grids shows that there 
is agreement at 2 and 4 times the original linear resolution. 

Vortices are visible in the inviscid runs for both planet masses 
ji = 10 -3 and 10 -4 in the grid codes, which induce a strong per- 
turbation to the tidal torque. The vortices in the upwind simulations 
have a larger amplitude and are more extended than in the shock- 
capturing results. The total torque acting on the planet excluding 
the material inside the Roche lobe agree in order of magnitude for 
Jupiter models. The torque results for Neptune have greater dis- 
persion, possibly due to incomplete clearing of the gap, but agree 
nevertheless in the final value within a factor 2. 

It has been observed that shock-capturing codes show a large 
amount of filamentary small-scale structure unseen in model re- 
sults obtained with other codes. This is especially true for both 
Direct-Eulerian implementations, AMRA and FLASH. In addi- 
tion, AMRA results show enhanced small-scale structure when 
compared to FLASH. Extensive comparison tests of the two im- 
plementations has shown that much of the observed differences is 
due to use of more selective dissipation algorithm in AMRA. (The 
so-called flattening algorithm in AMRA is based on Eqns. A.7- 
A.10 from IColella and Woodward! Jl984h while FLASH uses Eq. 
A.2). After adopting the simplified version of the flattening algo- 
rithm in AMRA, the results closely matched those obtained with 



the FLASH code. Adding a small amou nt of artificial viscosity 
with c oefficient of 0.1, as recommended bv lColella and W oodward 
resulted in a further reduction of filamentary structures and 
substantial reduction of the strength of vortices located at the gap 
edges. 

The upwind codes have a smooth disc structure and do not 
show filaments in the inviscid simulations. This may be due to the 
fact that shock-capturing codes have small intrinsic viscosity in our 
problem in cylindrical coordinat es, where flow is dominated by ad- 
vection in only one dimension. iBrvden et alJ ll999l) have shown 
that van Leer based codes in polar coordinates may have low in- 
trinsic viscosity comparable with shock-capturing methods. It has 
been checked that none of the above changes are needed in AMRA 
if the grid resolution is increased twice. In this case the solution is 
much smoother and the vortices at the gap edges decay faster. 

The Cartesian implementations produce results that are com- 
parable to the other codes but there are differences in the gap struc- 
ture due to the open inner boundary. The depleted density distribu- 
tion in the inner disc in FLASH-AP produces different torques but 
the torque contribution from the outer disc is consistent with the 
cylindrical grid codes. 

SPH codes predict the shape of the gap correctly but do not 
resolve well low density regions where the number of particles is 
small. In addition the spiral wakes are weaker, possibly due to SPH 
being more dissipative. The Balsara switch included in the SPH- 
TREE code is used to reduce the shear component of artificial vis- 
cosity but it may also smooth out the shocks. An advantage of SPH 
codes is that the geometry of the problem is well adapted to a La- 
grangian scheme and the algorithm implementation is simpler than 
for Eulerian codes. The planet can be treated as a regular parti- 
cle which accretes material. Furthermore, it is possible to follow 
the trajectory of individual fluid elements and study the accretion 
flows. SPH codes are computationally more expensive than Eule- 
rian codes at the same resolution. Our high resolution tests indicate 
that higher resolution is needed in the SPH simulations to obtain 
results comparable to the Eulerian grid codes. 

Possible future work includes the comparison of high reso- 
lution runs using multi-level meshes to investigate the gas flow 
close to the planet, the study of the orbital shift o f a free-moving 
plan et and 3 -dimen sional simulations (see e.g. iKlev et alJ l2001: 
D'Angelo et al. 2003). The convergence of the results with reso- 
lution needs to be studied in detail. 

In closing we would like to reiterate that computational work 
might be regarded as an experiment, rather than a simulation. We 
have shown that different codes can give slightly different results 
for the same physical problem. Reproducibility of experimental re- 
sults is fundamental to the scientific process, and this standard must 
be applied to those performed with computers. Before a computa- 
tional result can be regarded as reliable, it must be confirmed by an 
independent test with a different code. 
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APPENDIX A: ADVICE FOR OTHERS 

Putting together this comparison has been something of a "learn- 
ing experience" for all concerned. Although not strictly scientific, 
we would like to share our experiences with others who may be 
contemplating similar comparisons. 

As with many things, advance planning is the most impor- 
tant. So far is possible, decide in advance which quantities should 
be monitored, and how often this should be done. What should be 
checked every timestep (or so), and what is only required at much 
less frequent intervals? Storage requirements are relevant to this: 
for example, writing out the total mass in the simulation is a lot 
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cheaper (both in terms of space and time) than outputting the en- 
tire density field. Changing the output quantities at a later date will 
often involve re-running computations, which will delay matters. 

Careful attention should also be paid to the format of the sub- 
mitted files. Each code generally has its own output format. Those 
co-ordinating the comparison do not have time to pick through each 
of these - automated processing must be the goal. Make sure that 
the format is carefully specified (since if there are two mutually in- 
compatible ways of doing something, it is certain that results with 
both ways will be submitted). As an aside, for grid based results, 
it is probably more sensible to write out the indices of each cell, 
rather than the co-ordinates themselves: integers are exact. Supply 
the tables to convert indices to co-ordinates separately. 

Pay similar attention to the problem specification itself. Some 
flexibility will inevitably be needed, but try to keep this to a min- 
imum. Again, the authors' experience is that anything left vague 
will be done in different ways by different groups. 

Communication is also hugely important. In addition to set- 
ting up a mailing list, the authors were able to hold several short 
meetings, using funds provided by the EU. These were crucial to 
moving the project forward. Better still would have been to have 
held a longer workshop (perhaps a week) where everyone could 
gather, discuss and run their codes together. 

We hope that future groups will find our experiences useful in 
planning their own code comparisons. 
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